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ABSTRACT 

We study the evolution of the cold gas content of galaxies by splitting the interstellar medium 
into its atomic and molecular hydrogen components, using the galaxy formation model 
GALFORM in the ACDM framework. We calculate the molecular-to-atomic hydrogen mass 
ratio, H2/HI, in each galaxy using two different approaches; the pressure-based empirical re- 
lation of Blitz & Rosolowsky and the theoretical model of Krumholz, McKeee & Tumlinson, 
and apply them to consistently calculate the star formation rates of galaxies. We find that the 
model based on the Blitz & Rosolowsky law predicts an HI mass function, ^^C0(1 — 0) lu- 
minosity function, correlations between the H2/HI ratio and stellar and cold gas mass, and 
infrared-CO luminosity relation in good agreement with local and high redshift observations. 
The HI mass function evolves weakly with redshift, with the number density of high mass 
galaxies decreasing with increasing redshift. In the case of the H2 mass function, the number 
density of massive galaxies increases strongly from z — Q io z = 2, followed by weak evolu- 
tion up to z = 4. We also find that the H2 /HI ratio of galaxies is strongly dependent on stellar 
and cold gas mass, and also on redshift. The slopes of the correlations between H2/HI and 
stellar and cold gas mass hardly evolve, but the normalisation increases by up to two orders 
of magnitude from z = to z = 8. The strong evolution in the H2 mass function and the 
H2/HI ratio is primarily due to the evolution in the sizes of galaxies and secondarily, in the 
gas fractions. The predicted cosmic density evolution of HI agrees with the observed evolu- 
tion inferred from damped-Lya systems, and is always dominated by the HI content of low 
and intermediate mass halos. We find that previous theoretical studies have largely overesti- 
mated the redshift evolution of the global H2 /HI ratio due to limited resolution. We predict a 
maximum of pn^l Ptii ~ 1-2 at z « 3.5. 
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1 INTRODUCTION 

Star foimation (SF) is a key process in galaxy formation and evo- 
lution. A proper understanding of SF and the mechanisms reg- 
ulating it are necessary to reliably predict galaxy evolution. In 
recent years there has been a growing interest in modeling SF 
with sub-grid physics in cosmological scenarios, in which an 
accurate description of the interstellar m edium (ISM) of galax- 
ies is neede d (e.^ . ISpringel & Hernqui st 2003; Schave 2004; 



Dutton et aL* '2OIO'; 'Narayanan et al.' '2009; Schav e et all i2010l : 
Cook et al. 2010; Lagos et al. 201 1; Fu et al. 2010). 

It has been sh own observationally that SF ta kes place in 
molecular clouds (see lSolomon & Vanden Boutl2005l for a review). 
Moreover, the surface density of the star formation rate (SFR) 
correlates with the surface density of molecular hydr ogen, H2, in 
an approximately linear fa shion, Esfr oc Ehj (e.g. iBigiel et al.l 
I2OO8I : ISchruba et Zl l2010h . On the other hand, the correlation 
between Esfr and the surface density of atomic hydrogen, HI, 



is much weaker. Low SFRs have been measured in low stellar 
mass, Hl-dominated dwarf and low surface br ightness galaxies (e.g. 
iRovchowdhurv et alj|2009l : IWvder et a"Lll2009l) , and larger SFRs in 
more massive galaxies with mor e abundant H2, such as normal 
spiral and starburst galaxies (e.g. lKennicutlll998l : IWong & Blitzl 
2002). 

Observational constraints on the HI and H2 content of galax- 
ies are now available for increasingly large samples. For HI, accu- 
rate measurements of the 21 cm einission in large surveys of local 
galaxies have been presented by Zwaan et al. (2005) using the HI 
Parkes A ll-Sky Survey (HIPA SS; .Mever et al...2004.) and more re- 
cently bv lMartin et al. 1 2010|) using the Are cibo Legacy Fast ALFA 
Survey (ALFALFA: ' Giovanelli et al.ll2005h . From these surveys it 
has been possible to probe the HI mass function (MF) down to HI 
masses of A^hi ~ 10® and to estimate the global HI mass den- 
sity at 2; = 0, S7hi = 3.6 — 4.2 x 10~* in units of the present day 
critical density. These HI selected galaxies are also characterised 
by weaker clustering (e.g. lBasilakos et al.l2007tlMever et al.l2007h 
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than optically selected samples (e.g. iNorberg et al.ll20o"ll) . indicat- 
ing that local HI selected galaxies are preferentially found in lower 
mass haloes than their optical counterparts. At high redshift in- 
formation about HI is very limited since it is based mostly on 
absorption-lin e measurements in the spectra of quasi-stellar ob- 
ject (QSO: e.g.lPeroux et_a l.'2003':' Prochaska et al.l20q5HRao et al ' 



200^ ; iGuimaraes et all 12009: Note rdaeme et al.l 120091: see iRauchI 



19981 for a review of this technique). These observations suggest 
very little evolution of JIhi up to z ~ 5. Intensity mapping of the 
21 cm emission line is one of the most promising techniques to 
estimate HI mass abundances at high redshifts. This technique has 
recently been applied to estimate global HI densities at interme- 
diate redshifts {z < 0.8), and has given estimates in agreement 
with the ones inferred from absorp tion-li ne measurements ( e.g. 
IVerheiien et al.ll2007l : lL^et al.ll2007l 2009: IChang et al.l2010l) . 

To study H2, it is generally necessary to use emission from 
other molecules as tracers, since H2 lacks a dipole moment, mak- 
ing emission from this molecule extremely weak and hard to de- 
tect in interstellar gas, which is typically cold. The most commonly 
used proxy for H2 is the ^^CO molecule (hereafter 'C O'), which is 
the se cond most abundant molecule in the Universe. iKeres et al] 
( I2OO3I) reported the first attempt to derive the local luminosity 
function (LF) of CO(l — 0) (the lowest energy transition of the 
CO molecule) from which they inferred the H2 MF and the local 
f^Hs = 1.1 ± 0.4 X 10~* assuming a constant CO(l - 0)- 
H2 conversion factor. It has not yet been possible to estimate the 
cosmic H2 abundance at high redshift. However, a few detections 
of H 2 absorption in the lines-of-sight to QSOs have been reported 
(e.g. l Noterdaeme et al.l20()8l : lTumlinson et al.l2010l : ISrianand et al.l 
I2OICI) . as well as CO dete ctions in a large number of luminous 



star-forming galaxies (e.g. Greve et al. I I2OO5I: iGeachetakllioogl : 
iDaddi et al.ll20ld : lTacconi et al.ll2010h . 

Measurements of the HI and H2 mass content, as well as other 
galaxy properties, are available in relatively large samples of local 
galaxies (running into a few hundreds), allowing the characterisa- 
tion of scaling relations between the cold gas and the stellar mass 
content. From these samples it has been possible to determine that 
the molecular-to-atomic gas ratio correlates with stellar mass, and 
that there is an anti-correlation between the Hl-to-stellar mass ratio 
and stellar mass (e.g . .Bgthwell et al. 2009; Catinella et al. 2010: 
ISaintonge et al.l I2OI ll) . However, these correlations exhibit large 
scatter and are either subject to biases in the construction of obser- 
vational samples, such as inhomogeneity in the selection criteria, 
or sample a very narrow range of galaxy properties. 

The observational constraints on HI and H2 at higher redshifts 
will improve dramatically over the next decade with the next gen- 
eration of radio and sub-millimeter telescopes such as the Aus- 
tralian SKA Pathfinder (ASKAPj Johnston et al. 2008), the Ka- 
roo AiTay Telescope ( Me erK AT: |B 00 th et al. 2009) and the Square 
Kilometre Array (SKA: ISchilizzi et al. ,2008) which aim to de- 
tect 21 cm emi ssion from HI, and the Ataca ma Large Millimeter 
Array (ALMA; IWootten & Thompson"2009') and the Large Mil- 
limeter Telescope (LMT; Hughes et al. 2010) which are designed 
to detect emission from molecules. Here we investigate the pre- 
dictions of galaxy formation models for the evolution of the HI 
and H2 gas content of galaxies, taki ng advantage of the devel- 
opment of rea listic SF models (e.g. Mac Low & Klessen 2004'; 



Pelupessv et al.. 2006; Blitz & Ros olowsky 2006; Krumholz e t id. 
200s : IPelupessv & Papadopouloj I2OO9I : see iMcKee & OstriE3 



galaxies at the present day. Using such a model, it is reasonable to 
extend the predictions to follow the evolution of the gas contents of 
galaxies towards high redshift. 

Until recently, the ISM of galaxies in semi-ana lytic mod- 
els was treated as a single star-forming phase (e.g. Cole et aU 
2000"; Springel et al. ' 200ll : ICattaneo et al.l2008l:lLagos et al.ll2008l : 



20071 for a review). In our approach, a successful model is one that. 



at the same time, reproduces the observed stellar masses, luminosi- 
ties, morphologies and the atomic and molecular gas contents of 



Somerville et al. 2008). The first attempts to predict the separate HI 
and H2 contents of galaxies in semi-analytic models postprocessed 
the output of sing le phase ISM treatments t o add this informat ion 
a posteriori (e.g. lObreschkow et al.ll2009al : IPower et al.ll2O10l) . It 
was only very recently that a proper fully self-consistent treatment 
of the ISM and S F in galaxies throughout the cosmological calcula - 
tion was made (Cook et al.l2010l : lFu et alj2010l : lLagos et alj201lh . 
We show in this work that this consistent treatment is necessary to 

make progress in understanding the gas contents of gala xies. 

W e use the semi-analytical model GALFORM jCole et al] 
I2OOOI) in a A cold dark matter (ACDM) co smolo gy with the im- 
proved treatment of SF implemented bv lLagos et all ( I2OI l[) . which 
explicitly splits the hydrogen content in the ISM of galaxies into 
HI and H2. Our aims are (i) to study whether the models are able 
to predict HI and H2 MFs in agreement with the observed ones at 
2 = 0, (ii) to follow the evolution of the MFs towards high-redshift, 

(iii) to compare with the observational results described above and 

(iv) to study scaling relations of H2/HI with galaxy properties. By 
doing so, it is possible to establish which physical processes are 
responsible for the evolution of HI and H2 in galaxies. 

This paper is organized as follows. In !j2]we describe the main 
characteristics of the GALFORM model. In !j3]we present local uni- 
verse scaling relations and compare with available observations. In 
glwe present the local HI and H2 MFs, and the Infrared (IR)-CO 
luminosity relation. We also investigate which galaxies dominate 
the HI and H2 MFs and predict their evolution up to z = 8, and 
study the HI mass density measured in stacked samples of galaxies. 
In ij5] we present predictions for the scaling relations of the H2/HI 
ratio with galaxy properties and analyse the mechanisms which un- 
derly these relations. [j6]presents the evolution of the cosmic densi- 
ties of HI and H2, and compares with observations. We also deter- 
mine which range of halo mass dominates these densities. Finally, 
we discuss our results and present our conclusions in [jT] 



2 MODELLING THE TWO-PHASE COLD GAS IN 
GALAXIES 

We study the evolution of the cold gas content of galaxies by 
splitting the ISM into atomic and molecular hydrogen compo- 
nents i n the GAL FORM semi-analytical model of galaxy for- 
mation jCole et alj|2bQ 0: Benson et al. 2003; Baugh et al. 20q3; 
iBower et al.ll2006l : lBenson & Bowej|2010l : iLagos etal.12011 ). 

The model takes into account the main physical processes that 
shape the formation and evolution of galaxies. These are: (i) the 
collapse and merging of dark matter (DM) halos, (ii) the shock- 
heating and radiative cooling of gas inside DM halos, leading to 
the formation of galactic disks, (iii) quiescent star formation in 
galaxy disks, (iv) feedback from supemovae (SNe), from active 
galactic nucleus (AGN) heating and from photo-ionization of the 
inter-galactic medium (IGM), (v) chemical enrichment of stars and 
gas, and (vi) galaxy mergers driven by dynamical friction within 
common DM halos which can trigger bursts of SF, and lea d to the 
formation of spheroi ds (for a review of these ingredients see lBaughl 
|2006I :Ib ensonll2010l) . Galaxy luminosities are computed from the 
predicted star formation and chemical enrichment histories using 
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a stellar population synthesis model. Dust extinction at different 
wavelengths is calculated self-consistently from the gas and metal 
contents of each galaxy and the predicted scale lengths of the 
disk and bulge com ponents using a ra diative transfer model (see 
ICole et alfcOOd and lLacev et alj2011bh . 

We base our study on the scheme of Lagos et al. (201 1 ; here- 
after Lll) in which parameter-free SF laws were implemented in 
GALFORM. In the next four subsections we briefly describe the dark 
matter (DM) merger trees, the GALFORM models considered, the 
main features introduced by Lll into GALFORM, and the impor- 
tance of the treatment of the ISM made in this work. 



2.1 Dark matter halo merger trees 

GALFORM requires th e formation histori es of DM halos to model 
galaxy formation (see ICole et al1l2000l) . To generate these histo- 
ries we use an im proved version of the Monte-Carlo scheme o f 
ICole et al.l ilOOd) . which was derived by [Parkinson etall ( l2008t) . 
The Parkinson et al. scheme is tuned to match merger trees ex- 
tracted from the Millennium simulation of lSpringel et al-lfcoOSl) . In 
this approach, realizations of the merger histories of halos are gen- 
erated over a range of halo masses. The range of masses simulated 
changes with redshift in order to follow a representative sample of 
halos, which cover a similar range of abundances at each epoch. 

By using Monte-Carlo generated merger histories, we can ex- 
tend the range of halo masses considered beyond the resolution 
limit of the Millennium simulation, in which the smallest resolved 
halo mass is ~ lO^'^/i'^M© at all redshifts. This is necessary to 
make an accurate census of the global cold gas density of the uni- 
verse which is dominated by low m ass galaxies in low mass haloes 
jKim et alj2o"Tll : |Power et alfcoid : Lll). 

We adopt a minimum halo mass of A/haio = 5 x W^h~^ Mq 
at z = to enable us to predict cold gas mass str uctures down 
to the current observed limits (i.e. A/hi ~ 10^ Aftn. iMartin et al.l 
I2OIOI) . At higher redshifts, this lower limit is scaled with redshift 
to roughly track the evolution of the break in the halo mass func- 
tion, so that we simulate objects with a comparable range of space 
densities at each redshift. This allows us to follow a representative 
sample of dark matter halos, ensu ring that we resolve the structures 
rich in cold gas at every redshift. IPower et"al1 fcOIOl) showed that 
using a fixed tree resolution, as imposed by A'^-body simulations, 
can lead to a substantial underestimate of the gas content of the 
universe at z > 3. 

The cosmological parameters are input parameters for the 
galaxy formation model and influence the parameter values adopted 
to describe the galaxy formation physics. The two models used in 
this paper have somewhat different cosmological parameters for 
historical reasons; these cannot be homogenized without revisit- 
ing the choi ce of the galaxy fo rmation parameters. The parame- 
ters used in iBaugh et alj bOOSh are a present-day matter density 
of Qm = 0.3, a cosmological constant Q,a ~ 0.7, a baryon den- 
sity of fibaryons = 0.04, a Hubblc coustaut h = 0.7 and a power 
spectr um normalization of erg — 0.93. In the case of iBower et al.l 
( l2006l) . fim = 0.25, ST^A = 0.75, Hbaryons = 0.045, h = 0.73 and 
as = 0.9. 



2.2 Galaxy formation models 

We us e as starting points thel Baugh et al.l ( l2005h and lBower et al.l 
( l2006l) models, hereafter referred to as Bau05 and Bow06 respec- 
tively. The most important differences between these models reside 



in the mechanism used to suppress SF in massive galaxies and the 
form of the stellar initial mass function (IMF) adopted in starbursts. 
The Bau05 model invokes superwinds driven by SNe, which expel 
gas from the halo with a mass ejection rate proportional to the SFR. 
The Bow06 model includes a treatment of the heating of halo gas by 
AGN, that suppresses gas cooling in haloes where the central black 
hole (BH) has an Eddington luminosity which exceeds the cooling 
luminosity by a specified factor. The Bau05 model assumes that 
the IMF i s top-heavy during starbursts (driven by galaxy mergers) 
and has a lKermicutl Il983l) form during quiescent SF in galactic 
disks. The Bow06 model assumes a univer sal Kermicutt IMF (see 

1 Almeida et aT]|200l lOonzalez et al.ll2009l and lLacev et"ani201Ibl 
for a detailed comparison between the two models). Both mod- 
els give good agreement with the observed bj- and K-band LFs 
at 2 = 0, but only the Bow06 model agrees with the K-band LF at 

2 > 0.5, regardless of the SF law used (see Lll). The insensitiv- 
ity of the LF of galaxies in optical and near-infrared (IR) bands to 
the choice of the SF law was already noted by LI 1 and interpreted 
as the result of an interplay between the different channels of SF 
activity in the model (i.e. burst and quiescent modes). 

An important modification made in this work with respect to 
the original assumptions of the Bau05 and Bow06 m odels resides 
in the p arameters for the reionisation model, following iLacev et al.l 
( l201Ibl) . It is assumed that no gas is allowed to cool in haloes with a 
circul ar velocity below Vcrit at redshifts below ZreiorJlBensonetal 



'2003^ Taking into account recent simulations bv lOkamoto et al 



(2008.) and observa tional constraints on the reionisation redshift 
( Spergel et al.l2003l) . we adopt Krit = 30 kms~^ and z reion = 10, 
in contrast with th e values adopted by iBaugh et al.l ( l2005h and 
iBower et"al] ( |2006|) (Vcrit = 50 — 60kms~^ and Zreion = 6). Even 
though this affects the gas content of low mass halos, the changes 
are not significant for the results shown in this paper 

2.3 The interstellar medium and star formation in galaxies 

We use the SF scheme implemented in GALFORM by Lll. Lll 
tested different SF laws in which the neutral hydrogen in the ISM 
is split into HI and H2. T wo of the most promising are: (i) the 
iBlitz & Rosolow skv ( 2006) empirical SF law and (ii) the theoret- 
ical SF law of [krumholz et al. ( 2009). We briefly describe these 
below. 

(i) The empirical SF law of Blitz & Rosolowsky is of the form. 



SsFR = VSF fmol Sg 



(1) 



are the surface densities of SFR and total 
cold gas mass, respectively, vsf is the inverse of the SF timescale 



for the molecular gas and 



^mol / Egas is the molecular to 



total gas mass surface density ratio. The molecular and total gas 
contents include the contribution from helium, while HI and H2 
only include hydrogen (which in total corresponds to a fraction 
Xu ~ 0.74 of the overall cold gas mass). The ratio fmoi depends on 
the internal hydrostatic pressure as Ehj /Shi ~ fmoi/(fmoi — 1) = 
(Poxt/Po)"- Here ke is Boltzman n's constant and the values of 
lysF = 5.25 ± 2.5 x 10"^" yr"^ dLerov et al.ll2008l : iBigiel et"ai] 
l20Ilh.a = 0.92 ± . 07 aii dlog(Pn/kRfcm-^K]) = 4.54 ±0.07 
dBlitz & Rosolowskv! |200^ are derived from fits to the observa- 
tional data . To calculate Pcxt, we use the approximation from 
lElmegreenlJl993t) , 



(2) 



where Sgas and E, are the total surface densities of the gas 
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and stars, respectively, and ag and cr* are their respective verti- 
cal velocity dispersions. We assume a constant gas velocity dis- 

i iersion, tjg = 10 kms~^, following recent observational results 
LeroY et al.l2008l) . In the case of cr*, we assume self-gravity of the 
stellar disk in the vertical direction, a* = VTrGhTsT, where h-t is 
the disk height (see Lll for details). We estimate hi, by assuming 
that it is proportional to the radial exponential scale length of the 
disk, as observed in local spiral galaxies, with Rcfi /h, ~ 7.3 ±1.2 
jKregel et alj|2002l) . Hereafter, we will refer to the version of the 
model in which the Blitz & Rosolowsky empirical SF law is ap- 
plied by appending BR to the mode l name. 

(ii) In the theoretical SF law of lKrumholz et al.l j2009l) . fmoi in 
Eq.[T]depends on the balance between the dissociation of molecules 
due to the far-UV interstellar radiation, and their formation on the 
surfaces of dust grains, and i^sf the inverse of the time required to 
convert all of the gas in a cloud into stars. In their prescription, fmoi 
depends on the metallicity of the gas and the surface density of a 
molecular cloud, Ecompicx. Kmmholz et al. relate Scompicx to the 
smooth gas surface density, Sgas, through a clumping factor that is 
set to c = 5 to reproduce the observed Esfr — Egas relation. In 
summary, the Kmmholz et al. SF law is, 

Esfr = !/SF(Egas) /mol Sgas, (3) 

where 

/j^ s -0.33 

l'SF(Egas) = !^SF X f j for Egas < Eo 

/y, s 0.33 

= '^SP X f for Egas > Eo (4) 

Here i/§f = 3.8 x 10"^° yr"^ and Eo = 85 Mq pc"^ Hereafter, 
we will denote the version of the model where the Kmmholz, Mc- 
Kee & Tumlinson theoretical SF law is applied by KMT. 

We split the total hydrogen mass component into its atomic 
and molecular forms and calculate the SFR based on the molecular 
gas mass. It is important to note that in the BR SF law the inverse of 
the SF timescale for the molecular gas, i^sf (Eq.lTJ, is constant, and 
therefore the SFR directly depends only on the molecular content, 
and only indirectly on the total cold gas content through the disk 
pressure. However, in the KMT SF law, vsy is a function of cold 
gas surface density (Eq. |4}, and therefore the SFR does not only 
depend on the H2 content, but also on the cold gas surface density 
in a non-linear fashion. 

For starbursts the situation is less clear. Observational un- 
certainties, such as the conversion factor between CO and H2 
in starbursts, and the intrinsic compactness of star-forming re- 
gions, have not allowed a good characterisation of the SF law (e.g. 
[KennicutI 1998; Genzel et al. 2010; Combes et al. 2011). Theoret- 
ically, it has been suggested that the SF law in starbursts is differ- 
ent from that in normal star-forming galaxies: the relation between 
Eh2 /Ehi and gas pressure is expected to be dramatically differ- 
ent i n environments of very high gas densities typical of starbursts 
jPel upessv et al. 2006, Pelupessv & Papadopoulos 2009), where 
the ISM is predicted to be always dominated by H2 independently 
of the gas pressure. For these reasons we choose to apply the BR 
and KMT laws only during quiescent SF (fuelled by cooled gas ac- 
cretion into galactic disks) and retain the original SF prescription 
for starbursts (see Cole et al. 2000 and LI 1 for details). In the latter, 
the SF timescale is proportional to the bulge dynamical timescale 
above a minimum floor value and involves the whole cold gas con- 
tent of the galaxy, S FR = Mcom/tsf (see lOranato et al.ll2006l 
an d lLacevetal.l2008l for details). Throughout this work we assume 




Figure 1. Surface density of HI and H2 (solid and dashed lines, respec- 
tively) as a function of radius in units of the half-mass radius for randomly 
chosen galaxies in the Bow06 model, applying the KMT (top panel) and 
BR (bottom panel) SF laws. 



that in starbursts, the cold gas content is fully molecular, fmoi = 1- 
Note that this is similar to assuming that the BR pressure-law holds 
in starbursts (except with a different vsf) given that large gas and 
stellar densities lead to fmoi ~ 1. 

2.3. 1 Radial profiles of atomic and molecular hydrogen 

In order to visualize the behaviour of the HI and H2 components 
of the ISM of galaxies predicted in the models, we have taken the 
output of the original Bow06 model and postprocessed it to calcu- 
late the HI and H2 surface density profiles based on the expressions 
above. Fig. [T] shows the surface density profiles of HI and H2 for 
randomly chosen model spiral galaxies on applying the KMT (top 
panel) or the BR (bottom panel) laws. The HI extends to much 
larger radii than the H2. This is a consequence of the total gas sur- 
face density dependence of the H2 /HI ratio in both SF laws. Note 
that the KMT SF law gives a steeper decrease in the radial profile 
of H2 compared to the BR SF law. The BR SF law depends on the 
surface density of gas and stars, through the disk pressure, while 
the KMT SF law depends on the surface density and metallicity of 
the gas. 

We have compared our predictions for the size of the HI and 
H2 disks with observations. The models predict a coiTelation be- 
tween the HI isophotal radius, /hi, defined at a HI surface density of 
IM0 pc~^, and the enclosed HI mass, MHi(fei), which is in very 
good agreement with observations, irrespective of the SF law ap- 
plied. Both models predict A/hi(/hi)/Mq « 2x 10''(ZHi/kpc)^'^, 
while the observed relation f or Ursa Major is Mhi(/hi) /M0 ~ 
1.8 X 10^(ZHi/kpc)^ **^ (e.g. IVerheiien & Sancisill200ll) . On the 
other hand, the median of the relation between the exponential 
scale length of the H2 and HI disks predicted by the model is 
7?H2 ~ 0.4 Rm, whi le the relatio n infer red b y combining the ob- 
servat ional results of iRegan et al.l ( l200lh and IVerheiien & Sancisil 
1 I2OOII) is i?H2 = (0.44 ± 0.12) Rm- The model agrees well with 
the observed relations indicating that, in general terms, the ISM of 
modelled galaxies looks realistic. 
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Figure 2. The HI mass function at z = for the original Bow06 model 
when a constant Mjnoi/Mcoid = 0.27 ratio is assumed (dotted line), 
when a variable conversion factor based on the BR pressure law is assumed 
(dotted line) and for the Bow06.BR model (solid line), where the BR SF 
law is consistently applied throug hout the calculation . Symbols show ob- 
servational results at 2 = from IZwaan et al ] j2005l) using HIPASS and 
iMartin et al.. (2Qld) using ALFALFA, as labelled. 



2.4 Consistent calculation or postprocessing? 

Previous attempts to calculate the separate HI and H2 contents 
of galaxies in a cosmological scenario have been made by post- 
processing the output of existing semi-analytic models using 
specifically the B R SF prescription (e.g. Obreschkow et al. 2009a; 
iPower et al ]l2010h . Here we show that a self-consistent calculation 
of the ISM of galaxies, in which the new SF laws are included in 
the model, is necessary throughout in order to explain the observed 
gas properties of galaxies. To do this we choose as an example the 
Bow06 model and the BR SF law. 

Fig. |2] shows the HI MF in the original Bow06 model when 
a fixed c onversion ratio of Mmoi/Mcoid ~ 0.27 is assumed (dot- 
ted line; iBaugh et alj|2004h . when a variable conversion factor is 
calculated based on the BR pressure law (dashed line) and in the 
Bow06.BR model (solid line), where the BR SF law is applied 
consistently throughout the calculation. Symbols show a compila- 
tion of observational data. The original Bow06 model with a fixed 
H2/HI conversion gives a poor match to the observed HI MF. Post- 
processing the output of this model to implement a variable H2/HI 
ratio gives a different prediction which still disagrees with the ob- 
servations. If we consistently apply the BR SF law throughout the 
whole galaxy formation model, there is a substantial change in 
the model prediction which is also in much better agreement with 
the observations. These differences are mainly due to the fact that 
the SFR in the prescription of Bow06 scales linearly with the to- 
tal cold gas mass (as in starbursts, see §2.3), while in the case of 
the BR SF law this dependence is non-linear. Note that a thresh- 
old in gas surface density below which galaxies a r e not allowed 
to form stars, as originally proposed by lKennicut3 ( 1 19891) and as 
assumed in several othe r semi-analytic models (e.g. ICroton et al.l 
I2OO6I ; pRcce et alJl201Cll) . produces a much less pronounced, but 
still present dip (see LI 1). Thus, the linearity of the SFR with cold 



gas mass is the main driver of the strong dip at low HI masses in 
the original model. Remarkably, our new modelling helps repro- 
duce the observed number density even down to the current limits, 
A/hi ~ 10'' H^'^Mq. 

The new SF scheme to model the ISM of galaxies represents 
a step forward in understanding the gas content of galaxies. In the 
rest of the paper we make use of the models where the parameter- 
free SF laws are applied throughout the full calculation. 



3 SCALING RELATIONS FOR THE ATOMIC AND 
MOLECULAR CONTENTS OF GALAXIES IN THE 
LOCAL UNIVERSE 

Here we present the model predictions for various scaling relations 
between H2 and HI and other galaxy properties and compare with 
observations. In §3.1 we show how the H2/HI ratio scales with stel- 
lar and cold gas mass. In §3.2 we present predictions for this ratio 
as a function of galaxy morphology, and in §3.3, we study depen- 
dence of the HI and H2 masses on stellar mass. 

3.1 The dependence of H2/HI on galaxy mass 

The prescriptions described in §2 that split the ISM into its atomic 
and molecular components enable the model to directly predict cor- 
relations between the H2/HI ratio and galaxy properties. Fig. |3] 
shows the H2/HI ratio as a function of stellar mass (top panels) 
and total cold gas mass (bottom panels) at 2 = for the Bow06 
(solid lines) and the Bau05 (dashed lines) models using the BR and 
the KMT SF laws. The mo del predictions are compared to local 
observational estimates from lLerov et al.l ( l2008h . The errorbars on 
the model show the 10 and 90 percentiles of the model distribution 
in different mass bins. In the case of the observations, indicative er- 
rorbars on the H2/HI estimate due to the uncertainty in the CO-H2 
conversion factor are shown at the top of the top-left panel (i.e. the 
difference in the value inferred for starbursts and the Milky- Way, 
~ 0.6dex; see §4.2). Note that in this and subsequent figures where 
we compare with observations, we plot masses in units of h~^MQ 
to match observational units. For other plots we use the simulation 
units, h~'^MQ. 

Leroy et al. estimated stellar masses from 3.6/im luminosi- 
ties, which were transformed into /iT-band luminosities using an 
empirical conversion . Stella r masses were then calculated from the 
relation of lBell etal.U2003b between the stellar mass-to-light ratio 
in the Tf— band and B — V coloufl Thus the error on the stellar 
mass might be as large as a factor of 1.5 (considering the dispersion 
of 0.15 dex in the B — V vs. A'-band relation). 

The Bow06 model with the BR SF law predicts H2/HI ratios 
which are in very good agreement with the observed ones. The 
Bow06.KMT model fails to match the observed H2/HI ratios at 
all stellar and cold gas masses. This is due to the sharp decline of 
the radial profile of H2 (see Fig. which results in lower global 
H2/HI ratios in disagreement with the observed ones. In the case of 
the Bow06.BR model, it might appear unsurprising that this model 
is in good agreement with the observed correlation given that it is 
based on the empirical correlation between H2/HI and the hydro- 
static pressure in the disk (see §2). However, it is only because we 

^ Given that Leroy et al. adopt a i Kroup j feOOl!) IMF, we apply a small 
c orrection of 1.02 to adapt their stellar masses to account for our choice of 
a lKennicutll il983h IMF 



6 Claudia del P. Lagos 




7 8 9 10 11 7 8 9 10 11 



^fjLeroy et al. (2008) 




7 8 9 10 8 9 10 

log(M,„,yh^M^ log(M„,yh^M^ 



Figure 3. Moleculai-to-atomic hydrogen ratio, H2/HI, as a function of stel- 
lar mass (top panels) and total cold gas mass (M^oid = He+HI + H2; bot- 
tom panels) for the two model variants, the Bow06 (solid lines) and Bau05 
(dashed lines) models, when applying the BR (left panels) and the KMT 
(right panel) SF laws. Lines show the medians of the predicted distribu- 
tions, while errorbars show the 10% and 90% percentiles. Grey symbols 
show observations of spiral and in'egular galaxies from Leroy et al. ( 200^. 
Symbols with aiTows represent upper limits in the observed sample. In- 
dicative errorbars in the observational sample are shown in the bottom-left 
panel. 



are able to reproduce other galaxy properties, such as stellar mass 
functions (see §2), gas fractions and, approximately galaxy sizes 
(Lll), that we also predict the observed trend in the scaling rela- 
tions shown in Fig.[3l 

Both the Bau05.BR and BauOS.KMT models also give good 
agreement with the observed H2/HI ratios. The success of the 
BauOS.KMT model, in contrast with the Bow06.KMT model, is 
mainly due to the higher gas masses and metallicites predicted in 
the former model. However, the Bau05.BR and Bau05.KMT mod- 
els fail to reproduce the evolution of the K-hwA LF and the z = Q 
gas-to-luminosity ratios (see Lll for a complete discussion of the 
impact of each SF law on the two models). 

3.2 The dependence of H2/HI on galaxy morphology 

It has been shown observationally that the ratio of H2/HI masses 
correlates strongly with morphological type, with early-type galax- 
ies characterised by higher H2/HI ratios than late-type galaxies (e.g . 
lYoung & Kneze3ll989l : iBettoni et alj|2003l ; iLisenfeld et al]|201lh . 
Fig.|4]shows the H2/HI mass ratio as a function of the bulge-to-total 
luminosity ratio in the B-band, B /T, for different B-band absolute 
magnitude ranges for the Bow06.BR model. The right hand panel 
of Fig. [4] which shows the brightest galaxies, compares the model 
predictions with observations. All observational data have morpho- 
logical types derive d from a visual classification of the _B-band 
images ( de Vauco uleurs et al.l Q9 91), and have also been selected 
in blu e bands (e.g. Isimien & de Vaucouleurslll986l : I Weinzirl et al.l 
l2009i) . The comparison with observational data is shown for galax- 
ies in the model with Mb — 5 log(/i) < —19, which roughly coiTe- 
sponds to the selection criteria applied in the observational data. 



Note that we have plotted only galaxies in the model that have 
Mhi/M* ^ 10"^ and Myi^/M^, ^ 10~^, which correspond to 
the lowest HI and H2 gas fractions detected in the observational 
data shown. The Bow06.BR model predicts a relation between the 
H2/HI ratio and B /T that is in good agreement with the obser- 
vations. Note that, for the last bin, B /T < 0.2, the model pre- 
dicts slightly higher median H2/HI ratios than the values inferred 
from observations. However, in all cases a constant CO(l — 0)-H2 
conversion factor was assumed in the observational sample to in- 
fer the H2 mass. This might not be a good approximation in the 
low-metallicity environments typical of irregular or late -type spi- 
rals, where the H2 mass might be underestimated if a conversion 
factor typical of normal spiral galaxies is applied (e.g. lBoselli et al.l 
I2OO2I : see §4.2.1 for a discussion). 

The higher H2/HI ratios in early-type galaxies can be under- 
stood in the context of the dependence of the H2/HI ratio on gas 
pressure built into the BR SF law. Even though gas fractions in 
early-type galaxies are in general lower than in late-type galaxies, 
they are also systematically more compact tha n a late-type coun- 
terpart of the same mass (see lLagos et al.l201 lb . resulting in higher 
gas pressure, and consequently a higher H2/HI ratio. Also note that 
there is a dependence on galaxy luminosity: faint galaxies typically 
have lower H2/HI ratios than their bright counterparts. This is due 
to the contribution of the stellar surface density to the pressure, 
which increases in massive galaxies (see Eq|2j. 

3.3 The relation between HI, H2 and stellar mass 

Another form of scaling relation often studied observationally is 
the atomic or molecular hydrogen-to-stellar mass ratio as a func- 
tion of stellar mass. Recently, these relations have been reported 
for the atomic and molecular gas contents in a homogeneous sam- 
ple of relatively massive galaxies by Catinella et al. (2010) and 
Sain tonge et al.l ^201 ll) , respectively, with the aim of establishing 
fundamental relations between the stellar content of galaxies and 
their cold gas mass. Fig.|5]shows these relations for the Bow06.BR 
model z — compared to values reported for individual galaxies 
in these surveys with detected H2 or HI, respectivel50. Also shown 
in the top panel are a da ta compilation from the literature presented 
in iBothwell et l2009i) and individual galaxi es from the HERA 
CO Line Extragalactic Survey (HERACLES: iLerov et ai1l2009h . 
Measurements of AfH2 /A/stoiiar are subject to large errors (up to 
0.25 dex) given the uncertainty in the CO-H2 conversion factor (see 
discussion in §4.2.1). Observed Mni/Mstoiiar ratios are more ac- 
curate given the direct detection of HI. Horizontal lines in Fig. |5] 
show representative observational sensitivity limits (which are not 
exactly the same in all samples). Contours show the regions where 
different fractions of galaxies in the model lie, normalized in bins 
of stellar mass. The contours are not very sensitive to the location 
of the sensitivity limits. The model predicts the right location and 
scat ter of galaxies in thes e planes in contrast to previous models 
(see lSaintonge et alj201 ll for a discussion). Note that the literature 
compilation of Bothwell et al. shows larger scatter than the HER- 
ACLES, GALEX Arecibo SDSS Survey (GASS) and CO Legacy 
Database for the GASS survey (COLD GASS). Saintonge et al. 
suggest that this is due to the inhomogeneity of the literature com- 
pilation. 

Stellar masses in the observational samples were infen'cd using a 
Chabrier IMF. We scale them by a factor 0.89 to adapt them to our choice 
of a Kennicutt IMF. 
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Figure 4. Molecular-to-atomic hydrogen mass ratio, Mn^ / Mm, as a function of the bulge-to-total luminosity in the B-band, (B/T)b, in the Bow06.BR 
model, for galaxies with absolute B-band magnitudes, Mg — 5 log{h), in the range indicated on top of each panel. Contours show the regions within which 
different volume-weighted percentages of the galaxies lie for a given (B/T)g, with the scale shown by the key. For reference, the dashed line shows the 
median of the model distribution. Observational results from Young & Knezek 1 1989), Bettoni et al. 1 2003). Lerov et al. (2008) and Lisenfeld et aiHioT ih are 
shown as symbols in the right-hand panel, and we combine them so th at B/T < 0.2 corresponds to hi, Sm, Sd galaxies; 0.2 < B/T < 0.5 corresponds to 
Sc, Sb, Sa galaxies; B/T > 0.5 corresponds to E and SO galaxies (see lde Vaucouleurs et al.ll99ll for a description of each morphological type). 



The /A/stoiiar ratio in the model is only weakly corre- 
lated with stellar mass, in contrast to the A/Hi/Afstciiar ratio which 
is strongly dependent on stellar mass, in agreement with the ob- 
servations. The decrease in the MHi/Mstciiar ratio with increasing 
stellar mass is the dominant factor in determining the positive re- 
lation between H2/HI and stellar mass in Fig. |3] given the small 
variations of the /Msteiiar ratio with stellar mass. 

The predicted relation between and Af* (top-panel of 

Fig.[5j is close to linear for galaxies with Mh2 /Mi, > 0.05. These 
galaxies lie on the active star-forming sequence in the A/* — SFR 
plane (see Lll), due to an approximately linear relation between 
and SFR, although characterised by a large scatter (see 
§4.2.3). What drives the approximately constant SFR/Af* and 
A^H2/A^* ratios for galaxies on the active star-forming sequence 
is the balance between accretion and outflows, mainly regulated by 
the timescale for gas to be reincorporated into the host halo after 
ejection by SNe (see LI 1 for details). In the case of HI, the model 
predicts that, for galaxies plotted in the bottom-panel of Fig. [5] the 
HI weakly correlates with stellar mass, A/hi oc M^'^^, as a result 
of the feedback mechanisms included in the model. 



4 ATOMIC AND MOLECULAR HYDROGEN MASS 
FUNCTIONS 

The two-phase ISM scheme implemented in GALFORM (see §2) 
allows us to study the evolution of the HI and H2 in galaxies in 
terms of the MFs and their evolution. In the next two subsections 
we analyse the main mechanisms which shape the HI and H2 MFs 
and investigate how these interplay to determine the model predic- 
tions. 



4.1 Atomic hydrogen mass function 

Fig.|6]shows the z — OKI mass functions (HI MF) for the Bow06 
and Bau05 models u sing the BR and the KM T SF laws, and obse r- 
vational results from lZwaan et al.l(l2005h and lMartin et al.l ( l201oh . 



With the BR SF law, both the Bow06 and Bau05 models give 
predictions which are in reasonable agreement with the observed 
HI MF. However, when applying the KMT SF law both models 
give a poor match to the observed HI MF, with the Bau05 model 
underpredicting the number density of massive HI galaxies, whilst 
the Bow06 model greatly ovetpredicts the abundance of galaxies 
around A/hi ~ lO^/i~^M0. The latter is expected from the poor 
agreement between the H2/HI-stellar mass scaling relation pre- 
dicted by Bow06.KMT and the observed one (Fig.[3ll. In the case 
of BauOS.KMT, the poor agreement is expected from the low gas- 
to-luminosity ratios reported by LI 1. 

ICook et all ( |2010() showed that, by including the BR SF law 
self-consistently in their semi-analytic model, the number density 
of low HI mass galaxies is substantially increased, in agreement 
with our results. However, in their model this effect is still not 
enough to bring the predictions into agreement with the observed 
HI MF. The overproduction of low mass galaxies is also seen in 
their predicted optical LFs, suggesting that the difference resides 
in the treatment of SN feedb ack and reionisation (see Bow er et"ail 
l200dlBenson & Boweil2oT o). In the case of the KMT SF law, there 
are no published results using the full SF law in semi-analytical 
models. IFu et al. 1 1 I2OIOI) implemented the H2/HI metallicity- and 
gas surface density-dependence of the KMT SF law, but assumed a 
constant SF timescale for the molecular gas, instead of the original 
dependence on gas suiface density (see Eq.|4j. Furthermore, these 
authors compared their predictions to the observations of the HI 
MF only over the restricted range 10^ < Mm/h~'^MQ < W^°, 
one decade in mass compared to the five decades plotted in Fig.|6] 
so that it is hard to judge how well this KMT- like model really per- 
forms. 

Given that neither of the models using the KMT SF law predict 
the right HI MF or H2/HI scaling relations, and that the Bau05 
model fails to predict the right /("-band LF at 2; > 0.5, we now 
focus on the Bow06.BR model when presenting the predictions for 
the gas contents of galaxies at low and high redshifts. 
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Figure 5. Top panel: Molecular hydrogen-to-stellar mass ratio as a function 
of stellar mass for the Bow06.BR model at z = 0. The horizontal dotted 
line shows an approximate sensitivity limit below which CO(l-O) is not de- 
tected in the different molecular' surveys. Contours show the regions within 
which different volume-weighted percentages of the galaxies lie for a given 
stellar mass and above the sensitivity hmit, with the scale shown by the key. 
For reference, the dashed line shows the median of the model distribut ions. 
Observational data fr om the HERACLES s urvey jLerov et al1l20Q9l) . the 
C OLD GASS survey iSaintonge et ^1201 ih and the literature compilation 
of lBothwell et alJj2Q09t) are shown as symbols. Bottom panel: As in the top 
panel but for the atomic hydrogen -to-stellar ma ss ratio. Observational data 
from the GASS catalogue I Catinel la et al .l2010h and the literature compila- 
tion of Bothwell et al. are shown as symbols. The dotted line shows the HI 
sensitivity limit of the GASS catalogue. 



4.1.1 The composition of the HI mass function 

The low mass end of the z = HI MF {Mm < lO^/i'^Me), 
is dominated by satellite galaxies, as can be seen in the top panel 
of Fig.[71 Central galaxies contribute less to the number density in 
this mass range due to reionisation. Hot gas in halos with circu- 
lar velocity Kirc < 30kms~^ is not allowed to cool at 2 < 10, 
thereby suppressing the accretion of cold gas onto central galax- 
ies hosted by these hal os. These hal oes have masses typically 
Mhaio < 10" h~^MQ jBensorJlzoTob . This suggests that the HI 
content of galaxies in isolation may lead to constraints on reionisa- 
tion (Kim et al. in prep.). The HI MF of satellite galaxies does not 
show the dip at low HI masses observed in the HI MF of central 
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Figure 6. The HI mass function at z = for the Bau()5.KMT (dotted line), 
Bau05.BR (dot-dashed line), Bow06.KMT (dashed line) and Bow06.BR 
(solid line) models . Symbols show the z = observational esti mates from 
IZwaan et al.l j2005l) using HIPASS and from lMartin et al.l i2010l) using AL- 
FALFA, as labelled. 



galaxies because the satellites were mainly formed before reionisa- 
tion. Galaxies of intermediate and high HI mass mainly correspond 
to central galaxies. The predominance in the HI MF of central and 
satellite galaxies in the high and low mass ends, respectively, is 
independent of the model adopted. 

The Bow06.BR model slightly overpredicts the number den- 
sity of galaxies in the mass range 3 x lO'^/i^^MQ < Mhi < 
3 X lO**h,^^M0. This is mainly due to the slightly larger radii of the 
model galaxies compared to observations (see Lll), which leads 
to lower pressure within the galactic disk and therefore to a slight 
overestimate of the atomic hydrogen content. 

The bottom panel of Fig. |7] shows the contribution to the HI 
MF from galaxies hosted by halos of different masses. The HI MF 
at intermediate and high HI masses, i.e. Mm > 5 X 10^/i~^MQ, 
is dominated by galaxies hosted by low and intermediate mass DM 
halos. A/halo ;$ 10^^/i~^Mq, while lower HI mass galaxies are 
primarily satellites in higher mass halos. This scale in the DM halo 
mass (A/halo ~ lO^^/i^^M©) has been shown to be set by the ef- 
ficient suppression of SF in higher mass DM hal oes, mainly drive n 
by AGN feedback which shuts down gas cooling l lKim et alj201 ll ). 
In lower mass halos, in which AGN feedback does not suppress 
gas cooling, the cold gas content scales with the stellar mass of the 
galaxy and with the mass of the host DM halo. 



4.1.2 Evolution of the HI mass function 

Fig.[8]shows the evolution of the HI MF from z = 8 to 2 = 0. There 
is a high abundance of low HI mass galaxies at high redshift. This 
reduces with declining redshift as gas is depleted mainly through 
quiescent SF and starbursts. The number density of galaxies with 
low HI masses. A/hi ^ lO^/i^^M0, increases with redshift, being 
a factor 4 larger at 2 = 6 than at 2; = 0. The number density of HI 
galaxies with masses 10^/i~^Mq < Mhi ^ lO**/i~^M0 increases 
by a factor of 2 from 2 = to 2 = 1, with little evolution to 2 = 4. 
At 2 > 4, the number density of galaxies in this mass range drops 



Oh H 




6 7 8 9 10 11 

log(MH,/h" 

Oh H 




6 7 8 9 10 11 

log(MH|/h" Mg) 



Figure 7. The HI mass function at 2: = in tlie Bow06.BR model, dis- 
tinguishing between the contribution from central (solid line) and satellite 
(dashed line) galaxies in the top panel, and from galaxies hosted by DM 
hales of different masses, as labelled, in the bottom panel. 



again. The high-mass end of the HI MF, Mm > 5 x lO^/i'^M©, 
grows hierarchically until z ~ 1, increasing in number density by 
more than two orders of magnitude. The evolution of the high-mass 
end tracks the formation of more massive haloes in which gas can 
cool, until AGN heating becomes important (at z ~ 1). At 2 < 1 
the high-mass end does not show appreciable evolution. The HI 
MF from the break upwards in mass is dominated by galaxies in 
intermediate mass halos. These are less affected by AGN feedback, 
driving the hierarchical growth in the number density in the high 
HI mass range. Higher mass halos, i\fhaio ^ 10^^ H^^Mq, are 
subject to AGN feedback, which suppresses the cooling flow and 
reduces the cold gas reservoir in the central galaxies of these haloes. 
Lower mass haloes, A/haio ;$ 10^^ Mq, are susceptible to SNe 
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Figure 8. The HI mass functions for the Bow06.BR model at different red- 
shift s, as labelled. For reference, we show the z = observational results 
from lZwaan et al.l <2005.) and .Martin et al.. aOlO,) . 



feedback, which depletes the cold gas supply by heating the gas and 
returning it to the hot halo. 

Note that our choice of parameters for reionisation affects 
mainly the low mass end of the MF at high redshifts. At z — the 
abundance of galaxies with Mm ^ lO*fe~^M0 would be lower by 
a factor of 2 if a lower photoionisation redshift cut of Zrcion ~ 6 
was assumed. 

The characterisation of the HI MF at redshifts higher than 
z ~ will be a difficult task in future observations. The first 
measurements will come from stacking of stellar mass- or SFR- 
selec ted galaxy samples, as has been don e in the local Universe 
(e.g. Verh eiien et al.l2007l:lLah et al.ll200l 2009). Fig.|9]shows the 
cumulative HI mass per unit volume, pm, at different redshifts 
for samples of galaxies in the Bow06.BR model that have stellar 
masses or SFRs larger than Afstoiiar or SFR, respectively. The es- 
timated pm from source stacking can be directly compared to our 
predictions for the same lower stellar mass or SFR limit. The me- 
dian HI mass of galaxies as a function of stellar mass and SFR, is 
shown in the bottom-panels of Fig. |9] Errorbars represent the 10 
and 90 percentiles of the distributions, and are shown for three dif- 
ferent redshifts. The HI mass-stellar mass relation becomes shal- 
lower with increasing redshift, but the scatter around the median 
depends only weakly on redshift. The turnover of the median HI 
mass at Mstdiar ^ 10^h~^ Mq is produced by AGN feedback that 
efficiently suppresses any further gas cooling onto massive galax- 
ies and, consequently, their cold gas content is reduced. The HI 
mass is weakly correlated with SFR, particularly at high redshift. 
The median HI mass of galaxies spanning SFRs in the range plot- 
ted decreases with increasing redshift. This suggests that in order 
to detect low HI masses in observations at 2 « 0.5 — 2, stellar 
mass selected samples should be more effective than SFR selected 
samples. However, it would still be necessary to sample down to 
very low stellar masses (A/stdiar > 10* H'^Mq). Upcoming HI 
surveys using telescopes such as ASKAP and MeerKAT, will be 
able to probe down to these HI masses. 
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Figure 9. Top panels: cumulative HI mass density for galaxies in the 
Bow06.BR model with stellar masses (left panel) or SFRs (right panel) 
larger than a given value at different redshifts, as labelled. Bottom panels: 
median HI mass of galaxies as a function of stellar mass (left panel) and 
SFR (right panel). Errorbars correspond to the 10 and 90 percentiles of the 
distributions, and are shown for three different redshifts . 



4.2 Molecular hydrogen mass function 

The cold gas content is affected by SF, feedback processes, accre- 
tion of new cooled gas, but also by the evolution of galaxy sizes, 
given that our prescriptions to calculate the H2 abundance depend 
explicitly on the gas density (see 2.3). Our aim is to disentangle 
which processes primarily determine the evolution of H2 in galax- 



4.2.1 The present-day CO(l — 0) luminosity function 

Observationally, the most commonly used tracer of the H2 
molecule is the CO molecule, and in particular, the CO(l — 0) tran- 
sition which is emitted in the densest, coldest regions of the ISM, 
where the H2 is locked up. Given that in the model we estimate the 
H2 content, we use a conversion factor to estimate the C0(1 — 0) 
emission for a given abundance of H2, 



Jco/Kkms (5) 

Here A'^h2 is the column density of H2 and Jco is the integrated 
C0(1 — 0) line intensity per unit surface area. The value of X 
has been inferred observationally in a few galaxies, mainly through 
virial estimations. Typical estimates fo r normal spiral galaxies 
range between X 2.0 - 3.5 (e.g. lYoung & Scovillel Il99ll : 
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Figure 10. The CO(l - 0) luminosity function a t 2 = for t he Bow06.BR 
model compared to the observational estimates of l Keres etalj (,2003) for B- 
band (triangles) and a 60 [im (filled circles) selected samples of galaxies. To 
calculate the CO(l — 0) luminosities we assume a H2-to-CO conversion 
factor, X, where A^Ha/cm^^ = x X IQ-^O j^^/k km s^^. In the top 
panel, solid lines show the model predictions using fixed values of X, as 
labelled, and the dashed line shows the LF using the metallicity-dependent 
conversion factor, X(Z), of Boselh et al. (2002). In the bottom panel, the 
solid line shows the LF resulting from a conversion factor X = 3.5 and 
the different lines show the contribution of galaxies in different metallicity 
ranges, as labelled, for this choice of X. 



iBoselli et alJl2002l : iBlitz et alllioOTh . However, systematic varia- 
tions in the value of X are both, theoretically predicted and inferred 
observationally. For instance, theoretical calculations predict that 
low metallicities, characteristic of dwarf galaxies, or high densities 
and large optical depths of molecular clumps in starburst galaxies, 
should change X by a factor of up t o 5 — 10 in eith er direction 
(e.g. Bell et al. 2007; Meijerink et al. 2007|; lBavet eta l. 2009). Ob- 
servations of dwarf galaxies fav our larger conversion factors (e.g. 
X « 7; e.g. lBoseni et"ai]|2002h , wh ile the opposite holds i n star- 
burst galaxies (e.g. X « 0.5; e.g. iMeier & Tumeij|2004l) . This 
suggests a metallicity-dependent conversion factor X. However, 
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estimates of the coiTelation between X and metallicity in nearby 
galaxies vary significantly, from finding no correlation, when virial 
equilibrium of giant molecular clouds (GMCs) is assumed (e.g. 
I Young & Scovil le 1991 ; Blitz et al. 2007), to correlations as strong 
as X oc {Z/Zq)~^ , when the total gas content is inferred from the 
dust content on assumin g metallicity depe ndent dust-to-gas ratios 
(e.g. lGuelineT al. 1993; B oselli et alj|2002h . 

We estimate the CO(l — 0) LF by using different conversion 
factors favoured by the observational estimates described above. 
The top panel of Fig. 1 101 shows the C0(1 — 0) LF at 2: = when 
different constant conversion factors are assumed (i.e. independent 
of galaxy properties; solid lines). Observational estimates of the 
CO(l — 0) LF, made using a i?-band and a 60 ^m selected sample, 
are plotted using symbols jKeres et"n]|2003l) . The model slightly 
underestimates the number density at for X > 1, but gives good 
agreement at fainter and at brighter luminosities for sufficiently 
large values of X (such as the ones inferred in normal spiral galax- 
ies). In the predicted CO(l — 0) LF we include all galaxies with a 
Leo > 10^ Jy km/s Mpc'^, while the LFs from Keres et al. were 
inferred from samples of galaxies selected using 60 fim or B-band 
fluxes. These criteria might bias the LF towards galaxies with large 
amounts of dust or large recent SF. More data is needed from blind 
CO surveys in order to characterise the CO LF in non-biased sam- 
ples of galaxies. This will be possible with new instruments such 
as the LMT 

In order to illustrate how much our predictions for the CO ( 1 — 
0) LF at 2 = wary when adopting a metallicity-dependent 
conversion factor, X{Z), inferred independently from observa- 
tions, we also plot in the top panel of Fig. [TO] the LF when the 
X{Z) relation from lBoselU et all ( l2002h is adopted (dashed line), 
log(X) = 0.5toi - l.O2t;5:TOlog(Z/Z0). Note that this correla- 
tion was determined using a sample of 12 galaxies with CO(l — 0) 
luminosities in quite a narrow range, Leo ~ 5 x 10^ — 5 x 
10^ Jy km/s Mpc^. On adopting this conversion factor, the model 
largely underestimates the break of the LF. This is due to the contri- 
bution of galaxies with different metallicities to the C0(1 — 0) LF, 
as shown in the bottom panel of Fig.[TO] The faint-end is dominated 
by low-metallicity galaxies (Z < Zq/3), while high-metallicity 
galaxies (Z > Zq/S) dominate the bright-end. A smaller X 
for low-metallicity galaxies combined with a larger X for high- 
metallicity galaxies would give better agreement with the observed 
data. However, such a dependence of X on metallicity is opposite 
to that inferred by Boselli et al. 

By considering a dependence of X on metallicity alone we 
are ignoring possible variations with other physical properties 
which could influence the state of GMCs, such as the interstel- 
lar far-UV radiation field and variations in the c olumn density 
of gas (see for instance fPelupessv et al. 2006; Ba yet et al]|2009l ; 



Pelupessv & Papadopoulos 2009|; lPaDadoDoulosll20l^ Recently, 
Obreschkow et al. (,2009b) showed that by using a simple phe- 



nomenological model to calculate the luminosity of different CO 
transitions, which includes information about the ISM of galaxies, 
the C0(1 — 0) luminosity function of Fig. [TO] can be reproduced. 
However, this modelling introduces several extra free parameters 
into the model which, in most cases, are not well constrained by ob- 
servations. A more detailed calculation of the CO LF which takes 
into account the characteristics of the local ISM environment is be- 
yond the scope of this paper and is addressed in a forthcoming pa- 
per (Lagos et al. 2011, in prep.). 

For simplicity, in the next subsection we use a fixed C0(1 — 
0)-H2 conversion factor of X = 3.5 for galaxies undergoing qui- 
escent SF and X = 0.5 for those experiencing starbursts. In the 
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Figure 11. The H2 MF (top panel) and the CO(l - 0) LF (bottom panel) at 
different redshifts, as labelled, for the Bow06.BR model. Two fixed CO(l — 
0)-H2 conversion factors are assumed: X = 3.5 for quiescent SF and X = 
0.5 for burst SF. For reference, we show the CO(l — 0) LFs estimated 
at z = by .Keres et aL t2003i) from B-band and 60 /im-selected galaxy 
samples. 



case of galaxies undergoing both SF modes, we use X = 3.5 and 
X = 0.5 for the quiescent and the burst components, respectively, 
following the discussion above. 



4.2.2 Evolution of the H2 mass function 

At high-redshift, measurements of the C0(J — > J — 1) LF are not 
currently available. ALMA will, however, provide measurements 
of molecular emission lines in high redshift galaxies with high ac- 
curacy. We therefore present in Fig. [TT]predictions for the H2 MF 
and CO(l - 0) LF up to z = 8. 

The high-mass end of the H2 MF shows strong evolution from 
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z = 8 to 2; = 4, with the number density of galaxies increasing 
by an order of magnitude. In contrast, the number density of low 
H2 mass galaxies stays approximately constant over the same red- 
shift range. From 2 = 4 to 2 = 2 the H2 MF hardly evolves, with 
the number density of galaxies remaining the same over the whole 
mass range. From 2 = 2 to 2 = 0, the number density of massive 
galaxies decreases by an order of magnitude, while the low-mass 
end decreases only by a factor ~ 3. The peak in the number den- 
sity of massive H2 galaxies at 2 = 2 — 3 coincid es with the peak 
of the SF activity (see Llh lFanidakis et ^llOlOh . in which huge 
amounts of H2 are consumed forming new stars. The following de- 
crease in the number density at 2 < 2 overlaps with strong galactic 
size evolution, where galaxies at lower redshift are systematically 
larger than their high redshift counterparts, reducing the gas sur- 
face density and therefore, the H2 fraction. We return to this point 
in ^ The peak in the number density of high H2 mass galaxies at 
2 = 2 — 3 and the following decrease, contrasts with the mono- 
tonic increase in the number density of high HI mass galaxies with 
time, suggesting a strong evolution of the H2/HI global ratio with 
redshift. We come back to this point in §6. 

The evolution of the CO(l — 0) LF with redshift is shown in 
the bottom panel of Fig. [TT] Note that the z = Q LF is very sim- 
ilar to the one shown in the bottom panel of Fig. [TO] in which we 
assume a fixed CO-H2 conversion X = 3.5 for all galaxies. This 
is due to the low number of starburst events at 2 = 0. However, 
at higher redshifts, starbursts contribute more to the LF; the most 
luminous events at high redshift, in terms of CO luminosity, corre- 
spond to starb ursts (see jj4.2.3). 

Recently lOeach et al.l ( I2OIII ) compared the evolution of the 
observed molecular-to-(stellar plus molecular mass) ratio (which 
approximates to the gas-to-baryonic ratio if the HI fraction is 
small), /gas = A/nioi/(A/stciiar + Mmol), with predictions for the 
Bow06.BR model at 2 ^ 2, and found that the model gives a good 
match to the observed /gas evolution after applying the same ob- 
servational selection cuts. 



4.2.3 The IR-CO luminosity relation 

One way to study the global relation between the SFR and the 
molecular gas mass in galaxies, and hence to constrain the SF law, 
is through the relation between the total IR luminosity, Lm, and the 
C0(1 — 0) luminosity, Leo ■ We here define the total IR luminosity 
to be the integral over the rest-frame wavelength range 8— 1000/^m, 
which approximates the total luminosity emitted by interstellar 
dust, free from contamination by starlight. In dusty star-forming 
galaxies, most of the UV radiation from young stars is absorbed by 
dust, and this dominates the heating of the dust, with only a small 
contribution to the heating from older stars. Under these conditions, 
and if there is no significant heating of the dust by an AGN, we ex- 
pect that Lm should be approximately proportional to the SFR, 
with a proportionality factor that depends mainly on the IMF. On 
the other hand, as already discussed, the CO(l — 0) luminosity has 
been found observationally to trace the molecular gas mass in local 
galaxies, albeit with a proportionality factor that is different in star- 
bursts from quiescently star-forming galaxies. Observations sug- 
gest that sub-millimeter galaxies (SMGs) and QSOs at high redshift 
lie on a similar IR-CO luminosity relation to luminous IR galaxies 
(LIRGs) and ultra-l uminous IR galaxies (ULIRGs) o bserved in the 
local Universe (see ISolomon & Vanden BoutlllOOSi for a review). 
We investigate here whether our model predictions are consistent 
with these observational results. 

We show in Fig.[T2]the predicted IR-C0(1 — 0) luminosity re- 



lation, Lm — L'qq, for the Bow06.BR model at different redshifts, 
compared to observational data for different types of galaxies. The 
predicted CO(l — 0) luminosities for the model galaxies are cal- 
culated from their H2 masses as in §4.2.1, using conversion factors 
X — 3.5 for quiescent galaxies and X = 0.5 for bursts. To fa- 
cilitate the comparison with observations, we here express the CO 
luminosities L'qq in units of Kkms~^ pc^, which corresponds to 
expressing the CO line intensity as a brightness temperature. We 
predict the I R luminosities of the mod el galaxies using the m ethod 
described in 



L luminosities ot tne mod el galaxies using tne m etnod 
Lacev et al.l boilbh and iGonzalez et all JioTlh (see 



also^ Lacev et al. 1 201 lah ). which uses a physical model for the dust 



extinction at each wavelength to calculate the total amount of stel- 
lar radiation absorbed by dust in each galaxy, which is then equal 
to its total IR luminosity. The dust model assumes a two-phase in- 
terstellar medium, with star-forming clouds embedded in a diffuse 
medium. The total mass of dust is predicted by GALFORM self- 
consistently from the cold gas mass and metallicity, assuming a 
dust-to-gas ratio which is proportional to the gas metallicity, while 
the radius of the diffuse dust component is assumed equal to that 
of the star-forming component, whether a quiescent disk or a burst, 
and is also predicted by GALFORM. 

We show the model predictions in Fig. [T2] separately for qui- 
escent (left panel) and starburst galaxies (right panel), where qui- 
escent galaxies are defined as those whose total SFR is dominated 
by SF taking place in the galactic disk (i.e. SFRdisk > SFRburst). 
Solid lines show the median of the predicted IR-CO relation at dif- 
ferent redshifts, while errorbars represent the 10 and 90 percentiles 
of the distributions. To facilitate the comparison between quiescent 
and starburst galaxies, the typical IR luminosity at a CO luminosity 
of 10^ h-'^Kkms-^ pc^ is shown as dotted lines in both pan- 
els. We see that the separate Lir — L'qq relations for quiescent 
and starburst galaxies depend only slightly on redshift, but that the 
relation for starbursts is offset to higher IR luminosities than for 
quiescent galaxies. There are two contributions to this offset in the 
model. The first is the different SF laws assumed in starbursts and 
in galaxy disks. By itself, this results in roughly 40 times larger 
IR luminosities at a given H2 mass for starbursts. However, as al- 
ready described, we also assume a CO-to-jy2 conversion factor X 
which is 7 times smaller in starbursts, which causes an offset in the 
Lib, — L'qq relation in the opposite sense. The combination of these 
two effects results in a net offset of roughly a factor of 6 in Lm at 
a given Leo- 

A similar bimodality in the IR- CO luminosit y plane has 
also been inferred observationally by iGenzel et all ilOVf ) and 
ICombes et al.l ( |201 1|) . However, these results rely on inferring the 
CO(l-O) luminosities from higher CO transitions when the lowest 
transitions are not available, which could be sign ificantly uncertain, 
as we discuss below (see also llvison et alj20I ih . 

For comparison, we also plot in Fig. [12] a selection of obser- 
vational data. We plot data for local LIRGs and for UV/optically- 
selected star-forming galaxies at 2 ~ 1 — 2 in the left panel, to 
compare with the model predictions for quiescent galaxies, and for 
local ULIRGs, SMGs at 2 ~ 1 - 3, and QSOs at 2 ~ - 6 in the 
right panel, to compare with the model predictions for starbursts. 
We note that there are important uncertainties in the observational 
data plotted for high redshift (2 > 1) objects. For these, the IR lu- 
minosities are actually inferred from observations at a single wave- 
length (24/im or 850/im), using an assum e d shap e for the SED of 
dust emission. (In addition, the iRiechersI ( I2OIII) data on 2 > 1 
QSOs are actually for FIR, i.e. 40 — 120/xm, luminosities, rather 
than total IR luminosities.) Furthermore, the CO(l — 0) luminosi- 
ties for 2 > 1 objects are in most cases also not direct measure- 
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Figure 12. Infrared luminosity as a function of the CO(l — 0) luminosity 
for tlie Bow06.BR model at 2 = (red lines), 2 = 2 (yellow lines) and 
2 = 6 (blue lines) for quiescent (left panel) and starburst galaxies (right 
panel). Solid lines show the median, while errorbars show the 10 and 90 
percentiles of the distributions. To aid the comparison between the quies- 
cent and starburst galaxies, dotted lines show the predicted median IR lumi- 
nosity at 2 = 2 for a CO luminosity of 10^ Kkms ~^ pc^. We also 
show the following observational data: local LIRGs from lGao & Solomoj 
<2004) (triangles); local ULIRGs from Solomon et al. (1997) (squares); 
QSOs at 2 5 0.15 from IScoville et al.. ( 2003) , Evans et al. ( 2006) and 
[Bertram et al.l j2007l) (asterisks); sta r-for ming galaxies at 2 1.2 and 
2 2.3 from lTacconi et al l bOlOh andlGenzel et alj ilOld) (triangles); 
SMGs at 1 < 2 < 3 fromlGreveetal.1 j2005h . ISoIomon & Vanden BoutI 
i2005l ) and'T^cc oni et alj |200^ (diamonds); and QSOs at 1.4 < 2 < 6.4 
from Riechers ( 201 ih (crosses). Note that most of the observational data on 
the CO(l— 0) luminosity at high redshifts are infeiTed from the luminosities 
of higher CO transitions rather than being directly measured. 



QSOs follow the same Lir-L^q relation as starburst galaxies (e.g. 
lEvans et al .l2006l : lRiecherj201l[) . The model is thus able to explain 
the Liu-L'qq relation for all objects without needing to include any 
heating of dust by AGN. This also agrees with previous theoreti- 
cal predictions which concluded that only higher CO transitions 
are affected by the presen ce of an AGN (e.g. lMeiierink et al.l2007l ; 
lObreschkow et al.ll2009bh . 



5 EVOLUTION OF SCALING RELATIONS OF THE H2 
TO HI RATIO 

Our model predicts that the HI and H2 MPs are characterised by 
radically different evolution with redshift. This implies strong evo- 
lution of the H2/HI ratio. In this section we analyse scaling relations 
of the H2/HI ratio with galaxy properties and track the evolution of 
these relations towards high redshift with the aim of understanding 
what drives them. 

Fig. [H] shows the median H2/HI ratio as a function of stellar 
mass, cold gas mass and halo mass at different redshifts for the 
Bow06.BR model. Errorbars indicate the 10 and 90 percentiles of 
the model distribution in different mass bins. In the right panel, the 
predictions for satellite and central galaxies are shown separately as 
dashed and solid lines, respectively. For clarity, percentile ranges 
are only shown for central galaxies in this panel. In the case of 
satellites, the spread around the median is usually larger than it is 
for central galaxies. 

Fig. [T3] shows that the H2/HI ratio correlates strongly with 
stellar and cold gas mass in an approximately power-law fashion. 
The normalisation of the coiTelation between H2/HI and stellar 
mass (left panel of Fig.[T3} evolves by two to three orders of mag- 
nitude from z — 8I0 z = towards smaller values. The evolution 
is milder (only ~ 1.5 dex) if one focuses instead on the correlation 
with cold gas mass. Interestingly, the slope of the correlation be- 
tween the H2/HI ratio and stellar or cold gas mass hardly evolves. 
The approximate scalings of the H2/HI mass ratio against stellar 
and cold gas mass are; 



ments, but are instead inferred from measurements of higher CO 
transitions J — >■ J — 1 (usually 4 — > 3 or 3 — >■ 2). The conversion 
from L'coiJ — > J — 1) to L'qq{1 — 0) is usually done assuming 
that the brightness temperature of the CO line is independent of the 
transition J —> J — 1, as would be the case if the CO lines are 
emitted from an optically thick medium in thermal equilibrium at a 
single temperature (as appears to be the case in local spiral galax- 
ies). In this case, the luminosity L'gn is inde pendent of the transi- 
tion studied jSolomon & Vanden Bout|[2005h . However, recent ob- 
servations have shown different brightness temper atures for differ- 
ent CO transitions in so me high-redshift galaxies dOanielson et al.l 
l20ld ; llvisonetalJ|20T l|). As a result, there could be large errors in 
the CO(l — 0) luminosities of high-redshift galaxies when they are 
inferred from higher transitions. 

Comparing the model predictions for the IR-CO relation with 
the observational data, we see that the predictions for quiescent 
galaxies at 2: = are in broad agreement with the observations 
of nearby LIRGs, while at 2; = 2 the model predicts partially 
the location of UV/optically-selected star-forming galaxies. In the 
case of starburst galaxies, the predicted relation agrees with the ob- 
servations of low-redshift ULIRGs and high-redshift SMGs, and 
also with the observations of QSOs at both low and high redshift. 
The latter is consistent with the suggestion from observations that 



Note that the relation with i\fstciiar is only valid up to z = 4. At 
higher redshifts, the slope of the correlation becomes shallower. 
The relation with Afcoid, however, has the same slope up to 2 = 
8. The normalisation of the relation of H2 /HI with Mstdiar has 
a stronger dependence on redshift compared to the relation with 
Afcoid- However, note that the distribution around the median is 
quite broad, as can be seen from the percentile range plotted in 
Fig. [13] so these expressions have be taken just as an illustration of 
the evolution of the model predictions. 

The trend between the H2 /HI ratio of galaxies and host halo 
mass depends strongly on whether central or satellite galaxies are 
considered. In the case of central galaxies (solid lines in the right 
panel of Fig. I13t. there is a correlation between the H2/HI ratio 
and halo mass that becomes shallower with decreasing redshift in 
intermediate and high mass halos (A/halo > 5 X 10"/i^^Mo). 
This change in slope is mainly due to the fact that at lower redshift 
(2 < 1.5), AGN feedback strongly suppresses SF in central galax- 
ies hosted by halos in this mass range, so that the stellar mass-halo 
mass correlation also becomes shallower and with an increasing 
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Figure 13. Molecular-to-atomic hydrogen ratio as a function of stellar mass (left panel), cold gas mass (middle panel) and DM halo mass (right panel) at 
different redshifts, as labelled, for the Bow06.BR model. Lines show the medians of the distributions, while errorbars show the 10 to 90 percentile ranges. For 
clarity, these eiTorbars are plotted only for three different redshifts. In the right panel, central and satellite galaxies are shown separately as solid and dashed 
lines, respectively. Errorbars are only shown for central galaxies in this panel. 



scatter. At high redshift, the stellar mass-halo mass correlation is 
tighter and steeper. In the case of satellite galaxies, the H2 /HI ra- 
tio does not vary with host halo mass, but exhibits a characteristic 
value that depends on redshift. The lower the redshift, the lower 
the characteristic H2 /HI ratio for satellites. Nonetheless, the spread 
around the median, in the case of satellite galaxies, is very large, i.e. 
of 1.5 dex at 2: = 0. The lack of correlation between the H2 /HI ra- 
tio and halo mass and the large scatter are due to the dependence of 
the H2 /HI ratio on galaxy properties (such as stellar and cold gas 
mass) rather than on the host halo mass. 

In order to disentangle what causes the strong evolution of 
the H2 /HI ratio with redshift, we study the evolution of the galaxy 
properties that are directly involved in the calculation of the H2 /HI 
ratio. These are the disk size, cold gas mass and stellar mass, which 
together determine the hydrostatic pressure of the disk. Fig. [14] 
shows the H2/HI ratio, the half-mass radius (rso), the cold gas 
mass, the stellar mass and the midplane hydrostatic pressure of the 
disk at 150, as functions of the cold baryonic mass of the galaxy, 
Mbar = Afstciiar + A/cold- These predictions are for the Bow06.BR 
model. Note that we only plot late-type galaxies (selected as those 
with a bulge-to-total stellar mass ratio B /T < 0.5) as these domi- 
nate the cold gas density at any redshift. 

In general, galaxies display strong evolution in size, mov- 
ing towards larger radii at lower redshifts. However, this evolution 
takes place independently of galaxy morphology, implying that the 
galactic size evolution is caused by size evolution of the host DM 
halos (i.e. driven by mergers with other DM halos and accretion of 
dark matter onto halos). Note that the evolution in radius is about an 
order of magnitude from 2 = 6 to 2 = 0, which easily explains the 
evolution by two order of magnitudes in the H2/HI ratio (which is 
-P°xP r"^'*'*). By increasing the disk size, the gas density and 
pressure decrease and therefore the molecular fraction decreases. 
The evolution in cold gas mass is mild enough so as not to signifi- 
cantly affect the H2/HI ratio. Observationally, galaxies of the same 
rest-frame UV luminosit y seem to be a factor ~ 2 — 3 smaller a t 
2 « 6 than at 2 « 2 (e.g. lBouwens et al.ll2004l : IOesch et alj|2010h . 
in good agreement with the factor of ~ 2 — 3 evolution in size 
predicted by the Bow06.BR model (see iLacev et aljr2011bl for a 



detailed comparison of sizes predicted by the model with the ob- 
served galaxy population). 

Fig. [14] shows that the stellar mass of galaxies, at a given 
baryonic mass, increases with decreasing redshift and, therefore, 
the corresponding cold gas mass decreases. This drives the hy- 
drostatic pressure to be reduced even further. The midplane hy- 
drostatic pressure evaluated at rso, fcxt(r5o)/kB (bottom panel of 
Fig.ll4t. evolves by more than two orders of magnitude at a given 
A/bar over the redshift range plotted. Note that the typical values of 
fcxt(r5o)/kB predicted at 2 = are comparable to those reported 
in observations of nearby galaxi es, Pext/fes = 10 ^ — 10^ cm~"^K 
(e.g. Blitz & Rosolowskv 2006; Lerov et al.ll2008l) . At higher red- 
shifts, the values of Pcxt(r5o) also overlap with the range in which 
the EHa/Sni-Pcxt correlation has been constrained observation- 
ally at 2 = 0. 

The evolution in galactic size is therefore the main factor re- 
sponsible for the predicted evolution in the H2/HI ratio at fixed 
baryonic mass in the model, with a minor contribution from other 
properties which also contribute to determining this quantity (i.e. 
cold gas and stellar mass). 



6 COSMIC EVOLUTION OF THE ATOMIC AND 
MOLECULAR GAS DENSITIES 

The H2/HI ratio depends strongly on galaxy properties and redshift. 
In this section we present predictions for the evolution of the global 
density of HI and H2 . 

The top panel of Fig. \\5\ shows the evolution of the global 
comoving mean density of all forms of hydrogen (solid line), HI 
(dashed line) and H2 (dot-dashed line), in units of the critical den- 
sity at 2 = 0, p/pc,z=o- Observational estimates of the HI and H2 
mass density at different redshifts and using different techniques 
are shown using symbols (see references in §1). If the reported val- 
ues of Phi and pn, include the contribution from helium, then we 
subtract this when plotting the data. 

The model predicts local universe gas densities in good agree- 
ment with the observed ones, which is expected from the good 
agreement with the HI MF and C0(1 - 0) LF. At high redshift. 
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Figure 14. Molecular-to-atomic hydrogen ratio (top panel), half-mass ra- 
dius, rso (second panel), cold gas mass (third panel), stellar mass (fourth 
panel) and midplane hydrostatic pressure at rso (bottom panel), as func- 
tions of the total cold baryonic mass in the galaxy (i.e. A/bar = A/stollar + 
Afcold) at different redshifts (as labelled) for the Bow06.BR model. Solid 
lines show the medians of the distributions, while errorbars show the 10 and 
90 percentiles, and for clarity are plotted only for three different redshifts. 
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Figure 15. Top panel: Global density of all forms of neutral hydrogen 
gas (solid line), atomic hydrogen (dashed line) and molecular hydrogen 
(dotted-dashed line) in units of the critical density at z = 0, as a func- 
tion of redshift for the Bow06.BR model. Observational estimates of the HI 
mass density plotte d are fromj Zwaan et al. (2005) and Martin et a lj <2010l) 
from t h e HI MF. andlPeroux et al. (2003), Prochaska et al. (20Q5). „Rao et alj 
<2006h . lGuimaraes et alj j2009l) andlNoterdaeme et al. (2009|) from DLAs. 
Also shown is the local Universe estimate of the H2 mass density from 
iKeres et all i2003h using the CO(l - 0) LF. Bottom panel: Global atomic- 
to-total neutral hydrogen (dot-dashed line), molecular-to-total neutral hy- 
drogen (dashed line) and molecular-to-atomic hydrogen (solid hne) mass 
ratios, as functions of redshift. For comparison, the evolution of the global 
molec ular-to-atomic hydrogen ratio predicted by lObreschkow & Rawlingj 
<2009l) is shown as a dotted line. 



the model predicts pHi in good agreement with the observed den- 
sity of HI inferred from damped-Lya systems (DLAs). Note that 
in our model, by definition, the HI is attached to galaxies ( s ee jj2 ). 
However, recent simulations bv'paucher-Giguere & Kerea (1201 ll) . 
lAltav et alJ (I2OIO) and FumagaUi et al. (2011) suggest that DLAs 
might not coiTespond exclusively to HI gas attached to galaxies, 
but also to HI clump s formed during the rapid accretion of cold 
gas (e.g. lBinnevlll977l : lHaehnelt et aljl99^ . Thus, the comparison 
between our predictions of pni and the values inferred from obser- 
vations from DLAs has to be considered with caution in mind. 

At 2 = the HI mass density is much larger than that in H2, 
due to the low molecular fractions in extended, low pressure galac- 
tic disks. At 2 ~ 2, where galactic disks are characterised by larger 
pressure and therefore, higher molecular fractions, the density of 
H2 slightly exceeds that of HI. This H2 domination extends up to 
2 ~ 5, above which HI again becomes the principle form of neutral 
hydrogen. These transitions are clearly seen in the bottom panel of 
Fig. [15] where the evolution of the HI/H, H2/H and H2/HI global 
ratios are shown. Note that the peak of the H2/HI global ratio is 
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Figure 16. Global density of atomic hydrogen (top panel) and molecular hy- 
drogen (bottom panel) in units of the critical density at z = 0, as a function 
of redshift for the Bow06.BR model. The contributions from different DM 
halo mass ranges are shown by different lines as labelled. Observational 
estimates of the total HI and H2 mass densities are also shown. Symbols 
correspond to different observations as labelled in Fig. 1151 



at 2: ~ 3.5. The predicted evolution of tlie H2/HI global ratio dif- 
fers greatly from the theoretical predictions reported previously by 
Obreschkow & Rawlings (2009, dotted line in the bottom panel of 
Fig.[T5l( and lPovyeretall ( |2010|) . which were obtained from post- 
processed semi-analytic models (see §2.4). These authors reported 
a monotonic increase of the H2/HI ratio with redshift, even up to 
global ratios of H2/HI~ 10. This difference with our results is due 
in part to resolution effects, since these papers used DM halos ex- 
tracted from the Millennium simulation (see §2.1), and therefore 
only sample haloes with Mhaio ^ lO^"h~^M0 at all redshifts. 
Consequently, these calculations were not able to resolve the galax- 
ies that dominate the HI global density at high redshift, inferring 
an artificially high global H2/HI ratio. We can confirm this asser- 
tion if we fix the halo mass resolution in the Bow06.BR model at 
Mhaio = 10^° h~^MQ, whereupon a global ratio of H2/HI^ 7 is 



attained at z = 8. This is in better agreem ent with the value pre- 
dicted bv lObreschkow & Rawling3 ( |2009|) . but is still ~ 1.5 times 
lower. We interpret this difference as being due to the postprocess- 
ing applied by Obreschkow et al. rather than the self-consistent cal- 
culation adopted here. We have already showed that postprocessing 
can lead to answers which differ substantially from the self consis- 
tent approach (see Fig|2j. We can summarise the predicted evolu- 
tion of the H2/HI global ratio, PH2/PHI, as being characterised by 
three stages. 



PH2 / Pk 



0.13(1 + 2)^-^ for z<2 

f for 2<z<4: 



(8) 



0.45(1 

3.7 {1 + z) 



-0.7 



for 



> 4. 



Fig.[T6]shows the contribution to pm (top panel) and (bot- 
tom panel) from halos of different mass. The overall evolution of 
pm is always dominated by low and intermediate mass halos with 
^4aio < 10^^ h~^MQ. At 2; > 1, the contribution from halos 
with masses of 10" /i"^Mq < Mhaio < 10^^ h^^Mo quickly 
drops, and ai z > 4 the same happens with halos in the mass range 
IO^^/i-^Mq < Mhaio < 10" h-^Mg. At higher redshifts, low 
mass halos become the primary hosts of HI mass. In contrast, the 
evolution of is always dominated by intermediate and high 
mass halos with Mum > W^^ h~^MQ, supporting our interpre- 
tation of the evolution of the HI and H2 MPs . This sugge s ts tha t 
the weak clustering of HI galaxies reported bv lMever et al.l ( UoOTh 
and Basila kos et al.l j2007l) is mainly due to the fact that Hl-selected 
galaxies are preferentially found in low mass halos, while in more 
massive halos the hydrogen content of galaxies has a larger contri- 
bution from H2. 

Note that the evolution of pm and pu2 with redshift is weaker 
than the evolution of the SFR density, psFR- As LI 1 reported, /Osfr 
increases by a factor of ~ 15 from 2; = to 2 = 3 (the predicted 
peak of the SF activity). On the other hand, pH2 evolves by a factor 
of ~ 7 from 2 = to 2 = 3 (where the H2 density peaks), while 
pm evolves weakly, decreasing by a factor only of 1.5 from 2 = 
to 2 = 5. The difference in the redshift evolution between and 
pSFB. is due to the contribution from starbursts to the latter at 2 > 2, 
that drops quickly at lower redshifts. This leads to a larger decrease 
in the global psFR compared to pu2 ■ We remind the reader that 
we assume different SF laws for starbursts and quiescent SF, and 
therefore, by construction, different gas depletion timescales (see 
§2.3). 

In order to gain an impression of how much the predicted 
global density evolution can change with the modelling of SF, 
Fig.|17|shows the global density evolution of the neutral H, HI and 
H2 for the Bau05.BR, BauOS.KMT, Bow06.BR and Bow06.KMT 
models. The former three models predict very similar global den- 
sity evolution of HI at 2 < 4, while at higher redshifts the 
Bau05.KMT model predicts larger global densities mainly due to 
the metallicity dependence of the H2/HI ratio in the KMT SF law 
(see Eq.|4ll and the typically lower gas metallicities in high redshift 
galaxies. However, the Bow06.KMT model predicts a global den- 
sity of HI offset by a factor of 4 with respect to the observational 
estimates and the other models over the whole redshift range. This 
is caused by the large overprediction of the number density of inter- 
mediate HI mass galaxies by this model (Fig. (6)1. In the case of the 
global density of H2, all the models predict a similar redshift peak 
at 2 ~ 3 — 4, but with an offset in the maximum density attained. 
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Figure 17. Global density of all fornis of neutral H (top panel), HI (mid- 
dle panel) and H2 (bottom panel) in units of the critical density at z = 0, 
as a function of redshift for the Bow06.BR, Bow06.KMT, Bau05.BR and 
BauOS.KMT models, as labelled in the top panel. Symbols show observa- 
tional estimates, as labelled in Fig. 1151 



7 SUMMARY AND CONCLUSIONS 

We have presented predictions for the HI and H2 content of galaxies 
and their evolution with redshift in the context of galaxy formation 
in the ACDM framework. We use the extension of the GALFORM 
semi-analytic model implemented by Lagos et al. (2010, Lll), in 
which the neutral hydrogen in the ISM of galaxies is split into its 
atomic and molecular components. To do that we have used ei- 
ther the empirical correlation of Blitz & Roslowsky (2006, BR), 
between the H2/HI ratio and the hydrostatic pressure of the disk 
or the formalism of Krumholz et al. (2009, KMT) of molecule 
formation on dust grains balanced by UV dissociation. We as- 
sume that the SFR depends directly on the H2 content of galax- 



rameter, so as to isolate the consequences of changing the SF law 
on the model predictions. 

We study the evolution of HI and H2 predicted by two mod- 
els, those of Baugh et al. (2005, Bau05) and Bower et al. (2006, 
Bow06), after applying the BR or the KMT SF law prescriptions. 
When applying the BR SF law, we find that the Bau05 and the 
Bow06 models predi ct a present-day HI MF in broad agreemen t 
with observations bvl Zwaan et al.l ( l2005h and lMartin et al.1 Jioioh . 
and a relation between the H2/HI ratio and ste llar and cold gas mass 
in broad agreement with iLerov et all ( 1200 When applying the 
KMT SF law, the two models fail to predict the observed HI MF at 
z = 0. Moreover, Lll find that the Bau05 model fails to predict the 
observed iC-band LF at z > 0.5, regardless of the SF law assumed. 
Hence, we select as our best model that of Bow06 with the BR 
SF law applied (Bow06.BR). This model predicts local Universe 
estimates of the HI and H2 MFs, scaling relations of the H2/HI, 

/Mateiiar and Mm I Mstoiiar ratios with galaxy properties, and 
also the IR-CO luminosity relation in reasonable agreement with 
observations, but also gives good agreement with stellar masses, 
optical LFs, gas-to-luminosity ratios, among others (see Lll for a 
full description of the pr edictions of the GALF ORM model when 
changing the SF law and iFanidakis et al]|2010l for predictions of 
this model for the AGN population). We show that the new SF law 
applied is a key model component to reproduce the local HI MF at 
intermediate and low HI masses given that it reduces SF and, there- 
fore, SNe feedback, allowing the survival of larger gas reservoirs. 

Using the Bow06.BR model, we can predict the evolution of 
the HI and H2 MFs. We find that the evolution of HI is charac- 
terised by a monotonic increase in the number density of massive 
galaxies with decreasing redshift. For H2, the number density of 
massive galaxies increases from 2; = 8 to z = 4, followed by very 
mild evolution down to 2: = 2. At 2: < 2, the number density in 
this range drops quickly. These radical differences in the HI and 
H2 MF evolution are due to the strong evolution in the H2/HI ra- 
tio with redshift, due to the dependence on galaxy properties. We 
also present predictions for the HI mass-stellar mass and HI mass- 
SFR relations, and for the cumulative HI mass density of stacked 
samples of galaxies selected by stellar mass or SFR at different red- 
shifts, which is relevant to the first observations of HI expected at 
high redshift. 

At a given redshift, the H2/HI ratio strongly depends on the 
stellar and cold gas masses. The normalisation of the correlation 
evolves strongly with redshift, increasing by two orders of magni- 
tude from z = up to 2: = 8, in contrast with the very mild evo- 
lution in the slope of the correlations. We also find that the H2/HI 
ratio correlates with halo mass only for central galaxies, where the 
slope and normalization of the coiTelation depends on redshift. For 
satellite galaxies, no correlation with the host halo mass is found, 
which is a result of the weak correlation between satellite galaxy 
properties and their host halo mass. We find that the main mecha- 
nism driving the evolution of the amplitude of the correlation be- 
tween the H2/HI ratio and stellar and cold gas mass, and with halo 
mass in the case of central galaxies, is the size evolution of galactic 
disks. Lower-redshift galaxies are roughly an order of magnitude 
larger than 2 = 8 galaxies with the same baryonic mass (cold gas 
and stellar mass). 

Finally, we studied the cosmic density evolution of HI and H2, 
and find that the former is characterised by very mild e volution with 
redshift, in agreement with observations of DLAs (e.g. lPeroux et al.l 
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I2OO3I ; iNoterdaeme et al.ll2009l) , while the density of H2 increases 
by a factor of 7 from z « to z ~ 3. We also predict that H2 
slightly dominates over the HI content of the universe between z ~ 
2 — 4 (with a peak of PH2/PHI ~ 1-2 at z ~ 3.5), and that it is 
preferentially found in intermediate mass DM halos. On the other 
hand, we find that the HI is mainly contained in low mass halos. The 
global H2/HI density ratio evolution with redshift is characterised 
by a rapid increase from z = to 2 = 2 as PH2/PHI oc (1 + z)^ '^, 
followed by a slow increase at 2 < z < 4 as pu^/ pm oc (1 + 
z)"'®, to finally decrease slowly with increasing redshift at z > 4 
as Phq/phi oc (1 + z)~°'^. We conclude that p revious studies 
( lObreschkow & Rawlingil2009l; I Power et al. overestimated 
this evolution due in part to halo mass resolution effects, as the 
galaxies which dominate the HI density at z > 2 were not resolved 
in these calculations. 

The results of this paper suggest radically different evolution 
of the cosmic densities of HI and H2, as well as a strong evolu- 
tion in the H2/HI ratio in galaxies. The next generation of radio 
and submm telescopes (such as ASKAP, MeerKAT, SKA, LMT 
and ALMA) will reveal the neutral gas content of the Universe and 
galaxies up to very high-redshifts, and will be able to test the pre- 
dictions made here. 
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